Melting of aluminium clusters 



E. G. Noya and J. P. K. Doye 
University Chemical Laboratory, Lensfield Road, Cambridge CB2 1EW, United Kingdom 

F. Calvo 

Laboratoire de Physique Quantique, IRSAMC, Universite Paul Sabatier, 
118 Route de Narbonne, F31062 Toulouse Cedex, France 
(Dated: February 2, 2008) 

The melting of Al clusters in the size range 49 < N < 62 has been studied using two model 
interatomic potentials. The results for the two models are significantly different. The glue potential 
exhibits a smooth relatively featureless heat capacity curve for all sizes except for N = 54 and 
N — 55, sizes at which icosahedral structures are favoured over the polytetrahedral. Gupta heat 
capacity curves, instead, show a well-defined peak that is indicative of a first-order-like transition. 
The differences between the two models reflect the different ground-state structures, and neither 
potential is able to reproduce or explain the size dependence of the melting transition recently 
observed in experiments. 



I. INTRODUCTION 

Cluster thermodynamics has been an active field of re- 
search during the last decades. It is now accepted that 
clusters do not melt at a singular temperature, but, in- 
stead, there is a finite range of temperatures over which 
solid-like and liquid-like isomers coexisti Nevertheless, 
in the cases for which a well-defined peak is observed 
in the heat capacity curves, melting is generally consid- 
ered as the finite size analogue of a first-order phase 
transition^ Except for some particular casesjf^ clus- 
ters usually melt at a temperature lower than the bulk 
melting temperature, and this temperature generally de- 
creases with the cluster size. The appearance of premelt- 
ing effects, evidenced by a small peak in the heat capacity 
curve or by a change of slope in the caloric curve before 
melting, is also fairly common. 

Recently, the melting behaviour of positively charged 
aluminium clusters with 49 to 62 atoms has been 
measured. 5 Except for N = 49 and 59, the heat capac- 
ity curves show well-defined peaks between 450 and 650 
K. The melting temperature has an irregular variation 
with the cluster size, suffering a sharp drop at N = 56, 
which was attributed to a structural transition. A plot 
of the latent heat against temperature resulted in peaks 
at sizes N = 51, 57 and 61. Moreover, premelting effects 
were observed 100 K before melting for Al^ and Al^~ 2 . In 
order to try to understand the origin of the features ob- 
served experimentally, a theoretical study of the thermal 
properties would be valuable. 

The existing literature on Al clusters has mainly 
focused on the structural and electronic proper- 
ties. However, in spite of the experimental 6 * 7 - 
and theoretical 8 ! 9 . 10 ! 11 ! 12 ! 13 . 14 ! 15 ! 16 ! 17 ! 18 ! 19 ! 20 ! 21 ! 22 . 23 ! 24 . 25 
effort, a consensus has not yet been achieved about the 
structure of many of the experimentally observed magic 
numbers. Only for large sizes (N > 250) the structure 
of Al clusters seems to have been rationalized. Mar- 
tin et al. showed that, above this size, magic numbers 
were due to geometric shell closings associated with face- 



centered-cubic octahedra (fee), 7 and this hypothesis is 
consistent with results of kinetic simulationsii 4 * 1 ^* 1 ^ For 
smaller clusters, ab initio calculations have been per- 
formed at selected sizes p^ii* 1 ^ 3 - Only at N = 13 was 
agreement on the structure reached, which is thought to 
be icosahedral AiLiS^ The structures at N = 55 and 
147 remain unassigned, as different theoretical calcula- 
tions found very different structures. In particular, for 
N = 55, icosahedraljiii decahedral^ cub-octahedral^ and 
disordered 9 - structures were found most stable, depending 
on the method used. 

The lower computational expense of empir- 
ical potentials have allowed more extensive 
searchsi 19 . 20 . 21 . 22 ! 23 ' 24 ! 25 However, the results show 
a strong dependence on the model potential used. The 
Murrell-Mottram potential 2 ^ predicts a competition 
between fee and icosahedral structures in the size 
regime N = 2-55, with strong magic numbers for 
the 38-atom truncated octahedron and the 54-atom 
uncentered Mackay icosahedronji&iiL 2 !! Similar results 
were obtained with Voter-Chen 2 ! and Gupta 2 ^ models. 
Both show special stabilities also for the 38-atom trun- 
cated octahedron and for the Mackay icosahedron 24,25 , 
except that the Gupta potential favours the complete 
55-atom icosahedron^ Very different structures were 
found, however, with the Sutton-Chen 2 ^ and glue 30 
potentials. The former favours unusual, somewhat 
disordered structures, that are a hybrid of close packed 
structures, decahedra and Mackay icosahedra, 21,31,32 
A still different set of structures were found with the 
glue model, which favours polytetrahedral structures in 
the whole size range, except for those sizes close to the 
complete Mackay icosahedrai 2 ^ 

It is the aim of the present paper to use simulation 
techniques to study the melting of Al clusters in the size 
range studied experimentally,- in particular to see if the 
experimental trends can be captured by semi-empirical 
potentials. Besides, reliably finding the global minima 
and performing long enough simulations to achieve well- 
converged thermodynamic results is only tractable for 



such potentials. Several potentials that were fitted to Al 
properties have been suggestedi22i2L22i22i31 Of these, we 
have chosen to use the Gupta^S and glue^fi potentials, as 
they cover two very different set of structuresi^SiSiSi For 
sodium clusters, it has now been shown that the size de- 
pendence of the melting behaviour primarly reflects the 
geometric structure of the solid clusters^ 3 So, it is our 
hope that the behaviour exhibited by these two poten- 
tials will be representative of those that favour icosahe- 
dral and polytetrahedral structures, respectively, in this 
size range. It is noteworthy that high symmetry polyte- 
trahedral structures are possible (D^ at N = 51, D$h at 
N = 57, and T d at N = 61)21^ at the sizes for which the 
experimental latent heat showed a maximum, making it 
particularly interesting to study a potential that favours 
polytetrahedral structures. 

II. METHODS 

We performed canonical Monte Carlo (MC) simula- 
tions using the parallel tempering (PT) method 3 ^ to 
study the melting of Al clusters in the size range N = 
49 - 62. Our simulations consisted of 10 million MC 
steps, following an initial equilibration period of 1 million 
steps, for each of the 48 trajectories used, whose temper- 
atures ranged from 10 to 1000 K depending on the cluster 
size and on the model. All the trajectories were initial- 
ized with the ground-state structures for each size and 
modeli22i2i The lowest energy geometries for sizes not 
reported previously were obtained using a basin-hopping 
global optimization method. 37 Exchange among different 
temperatures was attempted with a probability of 10%. 
Evaporation or fragmentation of clusters was avoided by 
adding a spherical repulsive hard wall of radius tq + r, 
where r is the cluster ground-state radius (calculated 
as the distance from the center of mass to the furthest 
atom), and r the equilibrium bulk interatomic distance 
for Gupta and 3/2 of the dimer bond length for the glue 
model. These constraining radii were large enough not to 
affect the location of the main peak in the heat capacity 
curves. The data from these simulations was processed 
using the multihistogram technique^ 

The two potentials used belong to the family of 
embedded-atom potentials. They comprise a pair po- 
tential plus a many-body term that depends on the elec- 
tronic density at the atomic position. Therefore, the total 
energy can be written: 

i jyii i 

where 4>(rij) is a pair potential, nj is the interatomic 
distance between atoms i and j, F(p~i) is an embed- 
ding function, and pi is the electronic density at site i, 
which is usually approximated as the linear superposi- 
tion of the atomic densities of the rest of the atoms, i.e. 
Pi = p{ r ij)- The main difference between the two 

potentials is that Gupta assumes analytic forms for the 
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FIG. 1: Comparison of 4>{r), p(r) and F(p) for Al Gupta and 
glue potentials in the effective pair format. 



pair potential, the embedding function and the density, 28 
whereas, in the glue model, these functions are deter- 
mined by the fitting process^ The best way to compare 
the two potentials is in the effective pair format^ The 
forms for (j)(r) and F(p) are non-unique, and in the effec- 
tive pair format they are chosen so that the embedding 
function has a minimum at the value of p appropriate 
for the crystal. As can be seen from Fig. ^ there are 
considerable differences between the two potentials. For 
example, the pair potential for the glue potential is sig- 
nificantly shallower, and shows a double well structure as 
compared to the single well in the Gupta pair potential. 

The potentials also differ in the number of proper- 
ties and configurations that have been used in the fit- 
ting process. The glue model was adjusted to reproduce 
the forces obtained by first-principles calculations for a 
variety of environments, including surfaces, clusters, liq- 
uids and crystals^ However, the four free parameters 
in the Gupta potential have been simply fitted to some 
Al bulk properties, namely, the lattice parameters and 
elastic moduli. 39 * Therefore, the glue potential is more 
likely to work well in a variety of situations, and, for ex- 
ample, has been particularly successful in modelling the 
self-diffusion in Al4£ 

For both potentials some thermodynamic properties of 
the Al clusters have been studied previouslyfS 4 * 4 ^* 4 ^ 3 " but 
the only overlap with the size range considered here are 
at N = 54-56 for the Gupta potential. 24 i 41 i 42 
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FIG. 2: (Colour online) Ground state structures for Al (a) Gupta and (b) glue clusters with 49 to 62 atoms. For the 61-atom 
glue cluster, the disclination network is also depicted. 



III. RESULTS 

The structures of the global minima for this size range 
are depicted in FigGl For the Gupta potential, all 
the structures are based upon the Mackay icosahedron, 
whereas for the glue potential, they are all polytetra- 
hedral (i.e., the whole structure can be divided up into 
tetrahedra with atoms at the vertices), except at N — 
54 and 55, which are icosahedral. For the Gupta po- 
tential, the 55-atom Mackay icosahedron is most stable, 
whereas for the glue potential, the 55-atom and 61-atom 
structures are particularly stable. For polytetrahedral 
structures, it is generally preferred to have five tetra- 
hedra around a nearest-neighbour contact, but beyond 
a certain size, edges surrounded by six tetrahedra must 
also be present. The network formed by these sixfold 
edges, termed disclination network, provides an useful 
way to characterize polytetrahedral structures. For N = 
61 the tetrahedral disclination network is depicted in Fig. 
121 and most of the other polytetrahedral structures are 
based upon this structure. It is noteworthy that the pos- 
sible high-symmetry polytetrahedral structures at N = 
51 and 57, which involve a linear and a trigonal disclina- 
tion network, mentioned in the introduction as potential 
candidates for the experimental peaks in the latent heat, 



are not favoured by this potential. However, whether 
the 61-atom structure leads to a particularly large latent 
heat, should provide an indicator of whether this sugges- 
tion could be correct. 

That the two potentials show completely different 
structures does not necessarily mean that one or both 
are bad potentials, but illustrates how difficult it is for 
a potential to correctly predict a cluster's structure, be- 
cause to do so, the potential must be able to model a 
whole host of bulk and surface properties of the mate- 
rial correctly. Indeed, it is not uncommon for potentials 
that purport to model the same material to exhibit very 
different structures^ 

Figure shows the calculated canonical heat capaci- 
ties as a function of temperature for Al clusters with 49 
to 62 atoms using the Gupta and glue models. The heat 
capacity curves for the two potentials show few similari- 
ties. The Gupta heat capacities show a fairly well-defined 
peak for almost all sizes, as expected for clusters with 
icosahedral geometries, and the results for N = 54-56 
are consistent with previous results Interestingly, 
premelting effects are observed for the sizes N — 58-62. 

The glue potential, instead, predicts a smooth tran- 
sition from the solid to the liquid state without an as- 
sociated latent heat or peak in the heat capacity curve. 
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FIG. 3: Heat capacities as a function of temperature for 
aluminium clusters from 49 to 62 atoms, predicted by the 
glue (solid line) and Gupta (dashed line) models. 



The sizes N — 54 and 55, which have icosahedral ground- 
state structures, represent an exception, and their heat 
capacity curves are more similar to that for the Gupta 
potential, showing a well-defined peak albeit broader and 
at higher temperature than for Gupta. Similar behaviour 
was observed in a previous simulation study using the 
glue models Sun and Gong found a well-defined peak 




FIG. 4: Comparison of the experimental melting temper- 
atures (open circles) with those predicted by Gupta (closed 
circles) . 



in the heat capacity of the icosahedral clusters at N — 
13 and 147, but the polytetrahedral 43-atom cluster un- 
derwent a more continuous transition, without a clear 
peak in the heat capacity. Small peaks are found in the 
low temperature region of the heat capacity curves of 
the clusters AI53, AI56 and AI57, which are indicative of 
prcmclting effects. 

The fact that the two potentials predict significantly 
different melting behaviours is not surprising, given the 
differences in the structure of the lowest-energy clusters. 
The absence of a well-defined heat capacity peak for clus- 
ters with polytetrahedral global minima reflects the basic 
structural similarities between the ground-state and the 
molten clusters. It is well-known that simple liquids have 
substantial polytetrahedral character^ Hence, melting 
in these clusters is associated with the gradual occupa- 
tion of structurally similar isomers of higher and higher 
energy, rather than the more usual cooperative transi- 
tion between two sets of structures that have different 
energies and entropies. Similar behaviour is also seen for 
Lennard-Jones clusters with 25-30 atoms, as they have 
polytetrahedral global minimal A more cooperative, 
first-order-like transition only appears once the structure 
of the global minimum changes. Sun and Gong were right 
to describe such a continuous melting as more akin to the 
melting of a "glass" (or perhaps more properly an "ideal 
glass" as the clusters are always in equilibrium) than the 
melting of an ordered solid that has a fundamentally dif- 
ferent 'symmetry' to that of the liquid. 

The results obtained by either potential differ signifi- 
cantly from those obtained in the experiments of Breaux 
et al£ Of the two models studied, only the Gupta heat 
capacities show well-defined peaks for most sizes, as ob- 
served in experiments. The continuous nature of the 
melting transition for most of the glue Al clusters qual- 
itatively disagrees with experiments. Breaux et al. also 
reported a non-monotonic variation of the melting tem- 
perature with the cluster size, with a maximum at N — 
55 and a sharp drop for larger sizes. In Fig. the size 
dependence of the Gupta and experimental melting tem- 
peratures is compared. There is some agreement, but also 
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FIG. 5: A structural analysis of the melting of the (a) Gupta 
AI58 and (b) glue AI56 clusters. The percentage of quenches 
leading to the icosahedral ground-state structure (closed cir- 
cles), to structures which have absorbed some or all the extra 
atoms in the surface (diamonds), and to high energy isomers 
(open circles) is plotted. In (b) the structures considered are 
the polytetrahedral ground-state structure (closed circles), an 
icosahedral isomer (diamonds) and high energy isomers (open 
circles) . 



significant differences. As one would expect for clusters 
with Mackay-based structures, the Gupta clusters ex- 
hibit a maximum at the complete 55-Mackay icosahedron 
magic number, and the melting temperature decreases 
monotonically on either adding or removing atoms. How- 
ever, this maximum is broad and there is not the sharp 
decrease seen experimentally in going from N = 55 to 56. 
This broadness is probably due to compensating changes 
to both the energy and entropy of melting. 47 Instead, 
the latent heat of melting is a more structurally sensitive 
quantity^ The latent heats obtained from the Gupta 
clusters show a size dependence very similar to the one 
found for the melting temperature (not shown), and the 
experimentally observed peaks at N — 51, 57 and 61 
are not reproduced by our simulations. As for the sizes 
at which prcmclting effects occur, there is disagreement 
both between our results and the experiments, and be- 
tween the two models. In particular, the premelting ef- 
fects at N = 51 and 52 were not reproduced by either of 
the models. 

Even though experiments concern positively charged 
clusters and our calculations were performed for neutral 
clusters, it is unlikely that the discrepancies can be at- 
tributed to the charge. The effect that charge might have 
in such large clusters is probably small, as has been il- 



lustrated for Na clusters^ Rather, deficiencies in the 
potentials, and hence the structures they predict, are the 
more likely cause. 

In what follows we will analyse the origin of the pre- 
melting effects for both potentials. For that purpose, 
the microscopic behaviour of the system was followed by 
performing systematic quenches starting from each one 
of the trajectories. The AI58 Gupta and AI56 glue clus- 
ters were chosen as examples that show premelting effects 
with each potential. 

In Fig. [5] we have plotted the percentage of quenches 
that lead to a particular minimum for Gupta AI58. The 
ground-state structure of this cluster is a Mackay icosa- 
hedron, with the three adatoms located at anti-Mackay 
positions, all on the same face (Fig. 0). The percent- 
age of quenches that go to this structure is very high 
at low temperatures, but drops abruptly between 200 
and 300 K. Between approximately 250 K to 500 K, the 
lowest-energy minimum coexists with a family of isomers 
in which some or all the extra atoms, that lie above the 
surface in the ground state, are accommodated in the sur- 
face. Note that similar structures consisting of a Mackay 
icosahedron with one of the extra atoms absorbed in 
its surface are the ground-state structures for AI56 and 
AI57 (Fig. [2J. At temperatures higher than 550 K, high 
energy isomers dominate and the cluster is completely 
melted. Therefore, the first peak in the heat capacity of 
Gupta AI58 can be attributed to a surface reconstruction 
in which some or all the extra atoms accommodate in 
the surface. The same type of surface reconstruction was 
seen for the other Gupta clusters that showed premelt- 
ing effects, e.g. for AI57 an isomer in which the two extra 
atoms are accommodated in the surface becomes com- 
petitive above 200 K. Intuitively, these structures will 
become less favourable as the number of atoms in the 
third shell of the Mackay icosahedron increases, which is 
consistent with the fact that premelting features become 
more subtle for larger sizes. 

For all sizes at which the glue model exhibits premelt- 
ing effects, namely, N = 53, 56 and 57, the lowest-lying 
minima are polytetrahedral^ However, for AI56 an icosa- 
hedral isomer becomes competitive at 180 K, and approx- 
imately 50% of the quenches lead to this minimum up to 
approximately 580 K (see Fig. 0). Only above this tem- 
perature do high-energy isomers dominate. Similar anal- 
yses with analogous results were also performed for AI53 
and AI56. From the microscopic analysis it seems clear 
that the premelting effects observed in the glue model 
can be attributed to this solid-solid transition from poly- 
tetrahedral to icosahedral structures. 

Such a transition from polytetrahedral to icosahedral 
structures is somehow surprising. Even though icosahe- 
dral structures are energetically competitive for the near- 
est sizes to a shell closing^ the solid-solid transition can 
only be explained if the icosahedral isomers have higher 
vibrational entropy. However, for pair potentials at least, 
this is not the usual situation; instead, polytetrahedral 
structures usually have a higher vibrational entropy be- 
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cause of the greater internal strains i 4 * The situation is 
reversed here. For Al 56 , the geometric mean vibrational 
frequency is 1.0533 times larger for the polytetrahedral 
structure, which means that they have a lower vibrational 
entropy. An estimate of the solid-solid transition can 
hence be obtained within the harmonic approximation. 49 
Following this procedure, the solid-solid transition is pre- 
dicted to take place at 160 K, which is in very good agree- 
ment with the temperature of 168 K obtained from our 
PT simulations. 



IV. CONCLUSIONS 

We have studied the melting behaviour of Al clusters 
in the range N = 49-62 by means of PT MC simulations 
using two different models. The results are qualitatively 
and quantitatively different depending on the potential 
used, and neither of the potentials is able to reproduce 
the melting behaviour observed experimentally. 5 Many 
of the differences in the melting behaviour of the two 
models can be attributed to the different lowest-energy 
structures, which are mainly polytetrahedral for the glue 
potential and predominantly icosahedral for the Gupta 
potential. These results therefore suggest that neither 
polytetrahedral or icosahedral structures can explain the 
experimental size-dependence of the melting behaviour. 
Instead, the actual structure of Al clusters in this size 
range remains a mystery. Furthermore, we feel that none 
of the other semiempirical potentials are any more likely 



to be successful in reproducing the Al thermal proper- 
ties. An explanation of the experimental results would 
need, therefore, the application of models at a higher 
level of theory, but the computational demands of those 
techniques are significantly higher. 

Nevertheless, the results shown here have helped to un- 
derstand what kind of melting behaviour can be expected 
from potentials with different structural preferences. The 
melting transition for potentials that favour polytetra- 
hedral structures is more akin to an ideal glass transi- 
tion than the finite size equivalent of the bulk first-order 
melting transition. It is particularly interesting that this 
more continuous character of the melting transition even 
holds for sizes, such as N = 61, at which the polyte- 
trahedral global minimum is particularly stable. Also, 
due to the larger relative importance of the many-body 
term over the pair potential, unusual features, such as the 
solid-solid transition from polytetrahedral to icosahedral 
structure driven by the vibrational entropy, can be found. 
However, models that exhibit icosahedral structures will 
tend to exhibit well-defined heat capacity peaks akin to 
bulk melting. 
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